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Abstract 

We simulated the fourier transform of the correlation function of the 
Ising model in two and three dimensions using a single cluster algo- 
rithm with improved estimators. The simulations are in agreement 
with series expansion and the available exact results in d = 2, which 
shows, that the cluster algorithm can succesfully be applied for cor- 
relations. We show as a further result that our data do not support 
a hypothesis of Fisher that in any d = 2 lattice the fourier transform 
of the correlation function depends on the lattice generating function 
only. In d = 3 our simulation are again in agreement with the results 
from the series expansion, except for the amplitudes f±, where we find 
/+//_ = 2.06(1). 



1 Introduction 



A simple way to characterize the behaviour of spin variables of an Ising model 
consists in the two point correlation function g XjV . Except for the trivial d = 1 
case a general expression has not been found. In d = 2 are only the exact 
expressions for small and large separations known [[l] . Apart from that most 
information comes from high temperature expansions 0,0, |J. Universal 
quantities related to the critical behaviour with scaling dimension can be 
calculated from 4 field theory []5], J7|. Numerical simulations are difficult 
for several reasons. Direct simulation of g x y requires an computational effort 
increasing with N 2 , where N is the number of lattice points. Within periodic 
boundary conditions g is translational invariant and can be reduced to the 
fourier transform of go jX -Since this still needs an effort increasing with N, 
only lattices of rather modest linear extension can be treated for d > 1. 
For many applications the knowledge of the fourier transform along one or 
two directions is sufficient which renders us a feasible problem. Even within 
that restriction success depends on the algorithm. In general one may adopt 
a single spin update (Metropolis or Heatbath algorithm or a cluster 
algorithm || |10|| . For thermal quantities the smaller autocorrelation time of 
the latter is compensated by the better vectorization of the former algorithm 



ll|] . For the fourier transform of the correlation function one has to measure 
| J2x G zkx a x \ 2 . Since both the spin variable a x and the phase factor vary 
rapidly with x, the resulting fluctuations prevent any accurate measurement. 
The problem of fluctuation is particular striking approaching the critical 
temperature. The improved estimators in a cluster algorithm have much less 
variation and lead to more accurate results. The aim of the present paper 
is to show that the correlation function can be reliably determined with the 
single cluster algorithm of Wolff by comparing our results with known 
exact results in d = 2 and the results from series expansion in d = 2 and 



d = 3 0, EL 12 



The paper is organized in the following way. In Section |2] we collect the 
necessary formalism on the correlation function. The cluster algorithm is 
described in section |3] and section f| contains the results of our simulations. 
In section [5] we give our conclusions. 
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2 Correlation functions 



The connected pair correlation function g£ for an Ising model is defined as 

gt y =(<T x cr y )-m 2 e(T c -T) (1) 

x, y denote the sites of a d- dimensional simple cubic lattice of linear dimension 
L with periodic boundary conditions and (A) means the thermal average at 
temperature T taken with an energy E 

E = - 8\x-y\,lCT x (T y (2) 
where 6\ x \ i projects on next neighbours of x — 0. 

Below the critical temperature a magnetization term has to be subtracted 
from ((J x a y ) in order to obtain a reasonable behaviour at large distances 
\x — y\. g c xy can be considered as a L d ■ L d matrix. Usually its inverse (the 
proper vertex function in field theory) has a simpler behaviour. By means of 
translation invariance we can write 

(3) 

k is a d dimensional vector whose components ki are restricted to the first 
Billouin zone \ki\ < tt and have the form ^rij with integer n*. The positive 
function g(k) is simply the fourier transform of g^. For \k\ « 1 it has the 
expansion 



(l + ek 2 + ...) (4) 



9(k) xiT) 

where x(T) denotes the susceptibility in units of the Curie susceptibility and 

£(T) the effective correlation length. Near T ~ T c we expect the following 

jr. 



dependance on the scaling variable r = 1 — 1 



X(T) = C ± \t 

H(T) = f±\r\ 
m 



— ;/ 



Br". (5) 



The ± signs in the amplitudes refer to T < T c . In d = 2 the critical indices 
and the amplitudes C and B are known exactly |13|, |1^, [15], |T^|, but the 
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exact known f± refer to the true correlation length which is defined by 
the exponential decay of the correlation function lng XiV ~ — \x — y\/£, t for 
large separations and may differ especially for T < T c from £. Theoretical 
information on f± in d = 2 and d = 3 comes from series expansion [0, |3j and 
on the ratio /+//_ in d = 3 from renormalized perturbation theory ||. As 
discussed in the next section g(k) can be estimated by numerical simulations 
only at special values of the arguments. In the following we adopt the choices 
kn = ko(l, 1, ..1) (radial direction) and ki = &o(l,0, ..0) (linear direction). 
This corresponds to two functions gi(ko) and gn(k ) depending only on a 
scalar argument k^. 

gdh) = g(k L ) 

gn(ko) = g(k R ) (6) 

Due to the rotationally invariant form of eq(|J) x an d £ can be determined 
from one of the functions (^). In d = 3 we will restrict our numerical anal- 
ysis to the small k region. In d = 2 we can test in addition the following 
hypothesis stated by Fisher et. al. 0. Motivated by the results of the high 
temperature expansion they conjectured, that for all d = 2 Ising models g(k) 
is a function of the lattice generating function q(k) only which is given on a 
square lattice by 

q (k) = 2 (cos (h) + cos (k 2 )) . (7) 

Since g(q(k)) is a function of one variable only , g^ and g^ must satisfy the 
following relation 

9R (K) = 9l (k ) (8) 
provided k' is connected to k by 

sin 2 (ho/2) = 2sin 2 (k' /2) . (9) 

Even if eq(|J) would not be true, expansion of 1/g in powers of q consists in 
a very accurate representation of g c x y . Writing 

^ = l-Y.H n (T)(q(k)Y (10) 

9{ K ) n 

only H n with n < 2 are present up to o(l/T) 12 . Converting eq. (|10|) in position 
space we get from eq. @ 

(9%l = a o {T) S x>y — a\ (T) b\ x -y\,i 

a 2 (T) (Vy[,2 + 2V»|,VS) + ° (1/T)12 (U) 
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where 8\ x - y \ jr projects on points x, y having euclidean distance r. The high 
temperature expansion of ref [[| for a«(T) are given for convenience in ap- 
pendix A. 02 (T) is of order (^) 10 and therefore very small. If the expansion 
fllTD is valid, the coefficients a,- can be determined from low order moments 



of the function 



a (T) 
ai (T) 
a 2 (T) 



-E 
-E 

O T 

fed 



+ 



2 + cosA; , cosk 

gn(ko) 
cosk 



9l (fa) 



2L ^ ^ (fc ) 



1 ^ cos2k 
Lr 9L (fa) 



(12) 



Note, that the sum J2k is only one dimensional. For 1/T < 0.4 the first two 
terms in eq. flTTl) give a very accurate representation of (g c )~],, which may be 



useful in applications |L8 . 



3 Numerical methods 

A direct simulation of the disconnected correlation function 

9x,y = (&x<7y) (13) 

is prohibitive both from storage requirement and computational effort, even 
the latter can be reduced by Boolean operations (since a x are two valued 
quantities) and multi spin update Taking translation invariance within 
periodic boundary conditions into account only g 0ty is needed, or 

9o,y = ja E (°x> a y+x > ) (14) 

x' 

Having cured the storage problem the computational problem remains that 
eq.flnp requires 0(L d ) algebraic operations per single spin update. For the 
observables discussed in section |] not the full fourier transform of g is needed 
but rather gi(fa) or 9n(fa) with fa ^ 

^(*b) = p(lE^-*^-l 2 )- as) 
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The value at fco = has to be determined by an extrapolation k$ — > 0. The 
difference of this extrapolated value and the actual value of the r.h.s. of 
eq.(JH)|) gives m 2 L d . Since k in eq.(|T5]) can assume only L different values 
the computational effort can be tolerated. However, eq. ([15]) will still not lead 
to any meaningful result if one uses a conventional MC Method with single 
spin update. The reason are the cancellations in the sum involved in eq.flTB]) 
necessary to obtain a finite value for Qr^l- The corresponding fluctuations 
in each single measurement prevent any convergence of the average value in- 
side present computing facilities. These fluctuations can be avoided by using 
a cluster algorithm. Its merit is not so much the reduced autocorrelation 
time between different configurations because this is balanced by the better 
vectorization of the local algorithms [|11|| , but the existence of improved esti- 
mators for spin expectation values as (|T3|). For the single cluster method the 
estimator for the correlation function can be written as (see Appendix B) 



^(M = { 1 iEe-<^r). 

I s zee J 



(16) 



{ } means an average over all single clusters C and s denotes the number 
of sites in C. Comparing eq. fllB"!) and eq.(|l^) we notice the absence of the 
spin factor <j x which is due to the fact that all cluster spins have the same 
value. Moreover the sum over x extends only over the cluster instead of the 
whole lattice as in eq. (p~5j) . Since the size of a cluster in the high temperature 
phase is governed by the correlation length £, the phase factor will not vary 
substantialy for fco£ of 0(1) which is precisely the region we are interested in 
an application of eq. (|T6|) . Since the cancellations are only needed to render 
9r,l to be an decreasing function of k they are no problem at all. Note that 
this improvements hold only for cases where the magnetization is small. For 
T < T c and finite m a fraction 0(m) of the clusters will extend over the 
whole lattice and the same problems as in the single spin MC will reappear. 
Therefore we can use our method in the ordered phase in a narrow region of 
T c -T only. 

For an error estimate one has to keep in mind that the autocorrelation time 
t a for a cluster algorithm even much smaller than for single spin update 
may not be small absolutely. Therefore we divide the number of generated 
clusters into % blocks of size n s . In each block the average of any observable 
is taken. The final mean value of the observable and its error are calculated 
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from the different means of the blocks. This procedure will protect against 
autocorrelation times t a < n s . Typically we take a ratio n s /riB ~ 100 — 500. 
Also functions of observables (f.e.cj^ 1 ) including their error may be easily 
estimated. 



4 Results from numerical simulations 

In d = 2 exact expressions for the behaviour of g~ 1 {k) at small values of k 
are known. For a comparism we used the cluster algorithm on a 160 x 160 
lattice to simulate Typically we generated 10 5 clusters in the high 

temperature phase and 10 4 clusters in the ordered phase. The range of 
temperatures r = ±0.1, r = ±0.08, r = ±0.05 and r = ±0.03 was motivated 
because the correlation length must be small compared to the lattice size in 
order to avoid any finite size corrections. In Fig. 1 we show (A;)r~ 7 ' 4 as 
function of sin 2 {k/2) for various temperatures T > T c . Superimposed are 
straight line fits to points restricted to £ 2 /c 2 < 1. For r > 0.08 the data 
exhibit a linear behaviour in sin 2 k/2 in the whole k region. Approaching 
T c deviations from this behaviour become important in agreement with the 
high temperature expansion. gj}{k) exhibits the same behaviour. Since the 
statistical accuracy of g^ is much better and the extrapolation of gi to k — ► 
is compatible with the one obtained from g^ we discuss only the latter. The 
straight line fits yield us values for \ and £ 2 according eq.@. These fits 
are represented by the lines in Fig. 1. The resulting x an d £ are shown 
in Fig. 2 on a double log scale as function of r. The lines give the exact 
theoretical result |L5], [16], [17]], which agree well with the extrapolations of our 
data. The same analysis can be done for T < T c . Due to the non vanishing 
magnetization less statistic is available and the quality is poorer as compared 
to T > T c . The deviation from a straight line behaviour sets in much earlier 
than for T > T c , as can be seen from the big difference of data and the 
corresponding lines fitted to the small k region. From the fits we obtain \ 
and £ which are also shown in Fig. 2. Again we notice the agreement with the 
exact results 16], [17|]. The value of 5'^ 1 (0) can be translated via eq.(|l]) in a 



value for the magnetization. The values for m are divided by the theoretical 



value 13, 14 
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As one sees from Fig. 2, this ratio is compatible with 1 inside the small 
errors. In this case straight line fits to m are not very meaningful, since 
eq.flTTj) predicts a substantical deviation from the power law (|5|) in our range 
of t values. To demonstrate the statistical accuracy we fit our values of x 
and £ with the power laws (||). In fit A we fix the exponents to the theoretical 
values 7 = 7/4 and v — 1, in fit C we fit in addition the exponents 7 resp. v 
simultaneously to the data for T < T c . The resulting parameters are displayed 
in table 1 together with the known exact values, resp. the values from the 
series expansion for the effective correlation length f± from ref [fj, H|. Both 
agree within the statistics which shows that the cluster algorithm leads to 
reasonable values for the correlation function. 

In d = 2 we can test apart from the behaviour near T c also the results of 
the high temperature expansion. If the coefficients H n for n > 2 in the 
expansion of g~ l can be neglected , ao and ai can be determined from gi 
and gn. The experimental value of coefficients a and a\ are obtained by 
eq. (|T2|) . These values divided by the series expansions ( p4"D are shown in 
Fig. 4 as function of 1/T for various values of the lattice size L. Below 
1/T = 0.40 the first two coefficients in the expansion ([[1]) computed by the 
high temperature series give a good description of the correlation function 
which can be used in practical applications ||18|| . The deviations above 0.4 



are not due to finite size effects. They signal not so much a breakdown of 
the series expansion for a 0> i but rather the appearance of higher coefficients 
a n which spoil the relations (0) between g -1 and ao,i. This is corroborated 
by the values of 02 obtained by eq. (|T2]) which are zero below 1/T = 0.40 and 
quickly exceed the value expected from eq.(p4|) by an order of magnitude. 
Whereas finite size corrections are negligiable for leading effects discussed up 
to now, a test of the Fisher hypothesis g(k) being a function of the lattice 
generating function only is very sensitive on small finite size corrections. 
This is because any deviation from the relation ( |T2"D can be a small effect 
only and the cancellations necessary to suppress the low order 1/T terms in 
H n are no longer effective at order 1/L. We can control these corrections by 
the number of clusters with loop number |20| non zero (clusters wrapping 
around the lattice at least once). Requiring the fraction of those clusters 
to be less than 10 -6 we have to choose on a 160 x 160 lattice an inverse 
temperature not larger than 0.42. In order to make a small effect visible, we 
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divide both g L l (sin 2 ^) and g^i^sin 2 ^-) by a common linear factor R 

R = 0.005 + 2.63sm 2 V2 (18) 

obtained by a fit to gj} . In Fig. 5 gj}{sin 2 ^-) / R and g^i^sin 2,1 ^) / R are 
shown as function of sm 2 y. If g _1 would be linear in sin 2 ^, both data 
points would fall on a constant line at 1. The deviation from 1 signalizes the 
appearance of higher coefficients H n with n > 2 as expected from the series 
expansion. If the hypothesis ( |T2"| ) is true, both sets of data points have to 
agree, which is obviously not the case inside our errors. Since at the level 
of accuracy of 10 -3 finite size corrections can be neglected, our data are in 
disagreement with Fishers conjecture. 

For d = 3 we restrict ourselves to the behaviour near T c = 1/0.22165 | 2T| , ^2 



The computational effort increases rapidly with L, especially in the ordered 
phase at T < Tc. The data for a 40 3 lattice shown in Fig. 6 needed ~ 30h on 
a CRAY YMP. Since the correlation lengths are smaller in d — 3, we can cover 
the same range of r as in d = 2. As before the function guiko) = g(k Q , k , k ) 
has smaller errors than g^. We do not show a plot of gji, since its behaviour 
is very similar to the d = 2 case. Performing linear fits in the small k region 
we obtain x( r )>£( T ) an d m{r) shown in Fig. 6 on a double log scale. From 
the linear behaviour we see that the power laws predicted by eq.(|5|) are well 
satisfied. To obtain values for exponents and amplitudes we performed the 
following fits with or without correction to scaling: 

X ± = C ± \r\^(l+a ± \r\) (19) 
e± = f±\r\- u (l + a' ± \r\) (20) 
m = B (-rf (1 + b ± \r\) . (21) 

In principle the correction terms should be parametrized by nonlinear power 
laws, but within our statistics a determination of an exponent in the non- 
leading term is not possible. We checked that replacing |r| by |t| 1//2 in the 
correction does not change the values of exponents or amplitudes. In fit A 
we fix the exponents to the generally accepted values plj, [22], [L2|] and leave 
out the corrections which are allowed in fit B. The fits C and D are similar 
to A and B except that also the exponents 7,^ and (3 are free parameters. 
The resulting values and the corresponding x 2 /DoF are given in table II. In 
the fits to x an d tu correction terms are definiteley needed, whereas for £ 
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inclusion of the latter does not improve x 2 . In the case of the susceptibility 
variation of 7 does not reduce x 2 1 DoF. Therefore we believe that fit B is the 
most reliable one and is also shown in Fig. 6. The values for the amplitudes 
are in agreement with the series expansion (last column in table II). Variation 
of the exponent v for the correlation length leads to a value very close to the 
accepted value of v. Since including corrections in fit B did neither change 
the values of f± nor improve x 2 /DoF fit D should not be trusted so much. 
In this case we adopt fit C as the most reliable one. A common feature of 
fit A,B and C is a consistent 5% deviation in /_ from the series expansion. 
The universal ratio (from fit C) /+ / /_ 

/+//_ = 2.063(12) (22) 

disagrees with ref |12|, [7| , but is in agreement with a recent calculation using 
renormalized perturbation expansion in d = 3 || and with experimental data 
U . In the case of m fit D is definitely preferred due to the acchieved x 2 1 DoF . 
Both exponent (3 and the amplitude B are somewhat smaller as compared 
to the series expansion but are consistent within the errors. In table II our 
final values for the amplitudes are underlined. With the exception of /_ 
the cluster Monte Carlo method yielded values compatible with the series 
expansion with similar accuaracy. 

5 Conclusion 

Due to fluctuations the correlation function g in Ising models for d > 1 can- 
not be determined with local Monte Carlo methods except for rather modest 
lattice sizes. Since the improved estimators in a cluster algorithm exhibit 
much less fluctuation the fourier transform of g can be measured. We used 
the known results in d = 2 to demonstrate the validity of our method. As a 
byproduct we show that a hypothesis of Fisher may not be true. In d = 3 our 
simulations agree with the values expected from series expansion except for 
the amplitude /_ of the effective correlation length in the ordered phase. The 
universal ratio /+//- = 2.06(1) is in agreement with recent field theoretical 
estimates. 
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A High temperature expansion 



Taking the fourier transform of eq.([TT|) and comparing the coefficients of 



cos(kin), cos(k2m) with eq. flLOD one finds 

a = 1 - H (T) - AH 2 (T) 
ai = ffi(T) 

a 2 = H 2 (T). (23) 
Hi have been given in |2| in terms of powers of v — tgh(l/T) which leads to 



a Q = 1 + 4v 2 + 12t> 4 + Uv 6 + 188v 8 + 836v 10 + o{v 12 ) 
ai = v(l + v 2 + 5v 4 + 2lv 6 + 96v s + AOlv 10 ) + o(v 13 ) 
a 2 = Av 10 + o(v 12 ) (24) 

B The improved estimator 



The fourier transform of eq. (|14D is 



m = j- d Y.e lk(x - y) M (25) 

For a Swendsen Wang [fj cluster decomposition of the lattice we have 

(a x a y ) = A x (C)A y (C) (a c a C ') (26) 

c,c 

A X .(C) takes the value 1 if x £ C and otherwise. Noting that cluster spins 
are independent variables 

= J2A x (C)A y (C) (27) 

c 
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so 



^o = tUeie^t} 

L I c xeC ) 



(28) 



sw 



The 1/ L d J2c m the Swendsen Wang algorithm translates into 1/s for a single 

so 



cluster algorithm [03 
from which we obtain eq. (|16D . 



sec* 



(29) 
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table I 





A 


C 




7 


7/4 


1.753(4) 


7/4 


c + 


0.98(2) 


0.97(9) 


0.962582 




0.0247(3) 


0.025(3) 


0.025537 


V 


1 


0.99(3) 


1 


U 


0.558(4) 


0.55(2) 


0.56702(5)i 


f- 


0.16(3) 


0.18(2) 


0.175(5)6 



table II 





A 


B 


C 


D 


series tL^ 


7 


1.237 


1.237 


1.168(10) 


1.235(44) 


1.237 


c+ 


1.198(5) 


1.093(13) 


1.432(33) 


1.12(19) 


1.103(1) 




0.191(4) 


0.211(11) 


0.237(11) 


0.219(41) 


0.223(3) 


X 2 /DoF 


81/18 


18/16 


29/17 


18/15 




V 


0.629 


0.629 


0.628(1) 


0.618(6) 


0.629 


u 


0.4997(2) 


0.4995(8) 


0.501(2) 


0.522(11) 


0.496(4) 


f- 


0.2415(8) 


0.238(2) 


0.243(1) 


0.251(8) 


0.251(1) 


X 2 /DoF 


64/16 


62/14 


63/15 


59/13 




P 


0.330 


0.330 


0.301(2) 


0.319(5) 


0.330 


B 


1.632(2) 


1.686(4) 


1.478(9) 


1.608(37) 


1.71(2) 


X 2 /DoF 


263/11 


35/9 


35/10 


24/8 
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Figure captions 



Fig. 1: The inverse fourier transform of the correlation function gj}{k) 
above T c as function of sin 2 | for r = 0.1,0.08,0.05 and 0.03 (from bottom 
to top). The straight lines are fits to the data for £ 2 k 2 < 1. 

Fig. 2: d = 2: The susceptibility x? the effective correlation length £ and 
the magnetization m divided by eq.flTTD as function of \r\ on a double log 
scale. ± indicates the sign of T — T c . Solid lines represent the fit C (fit A) 
to £ (x) described in the text. The dotted line is the constant 1. 

Fig. 3: The inverse fourier transform of the correlation function gj}{k) 
below T c as function of sm 2 | for r = —0.1, —0.08, —0.05 and —0.03 (from 
bottom to top). The straight lines are fits to the first two k values with k ^ 0. 

Fig. 4: The coefficients a$ and a\ from eq.(|l|) divided by the series ex- 
pansion (^) for different values of the lattice size L as function of 1/T. The 
critical temperature is indicated by the arrow. 

Fig. 5: The linear <7l(^o) an d the radial gn{k' ) correlation function at 
1/T = 0.42 as function of sin 2 ^. If Fishers' hypothesis holds both have 
to be equal. 

Fig. 6: d = 3: The susceptibility x, the effective correlation length £ and the 
magnetisation m as function of |r| on a double log scale. ± indicate the sign 
of T — T c . The straight lines correspond to fit B (fit C, fit D) for x {£,, m )- 
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